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D ■ ABSTRACT 
tin 

In any coUisionless A^-body code, there is an optimal choice for the smoothing 
^ . parameter that minimizes the average error in the force evaluations. We show 

^ . how to compute the optimal softening length in a direct-summation code and 

demonstrate that it varies roughly as A^^^^'^. 
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1. Introduction 

The role of smoothing in an A^-body code depends on the nature of the system being 
modelled. When simulating a system containing a modest number of point masses, like an 



^ ! open star cluster, the goal is to reproduce the exact level of graininess that exists in the 
real cluster. Ideally no softening would be used, but a nonzero softening length is usually 
^ i included to eliminate the divergent forces that would result from a close encounter between 
^ [ two "stars." When modelling globular clusters, the softening length may be increased 
to reduce the graininess of the potential to the level expected in a system containing a 
larger number of particles than can easily be handled in the computer. When simulating a 
system containing finite-sized objects, such as a galaxy cluster, the softening length may be 
increased still more to approximate the physical size of a typical object. 

But A^-body codes are often used to study the evolution of systems, like galaxies, 
that contain extremely large numbers of particles. A code designed for such problems is 
sometimes called "coUisionless," since the system being modelled has a two-body relaxation 
time that greatly exceeds the elapsed time of the simulation. In a coUisionless A^-body code, 
the particles do not represent real objects; they are simply Monte-Carlo realizations of the 
underlying, smoother mass distribution. The purpose of the smoothing in a coUisionless 
code is to generate from the A^ particles an accurate estimate, at each time step, of the 
gravitational forces corresponding to the much smoother system being modelled. 
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Smoothing parameters in A'"-body codes are often chosen in an ad hoc way. Here 
we show that the softening length in a colhsionless code can be chosen objectively, so 
as to minimize the average errors in the force determinations. Too small a value for the 
smoothing parameter yields an estimate that is overly noisy, reflecting flnite-iV fluctuations 
in the forces. Too large a smoothing parameter reduces the noise but increases the error 
from the "bias," i.e. the systematic misrepresentation of the force due to the failure to 
resolve real features with scale lengths less than the softening length. In general, a unique 
value for the smoothing parameter may be found such that the combined error from these 
two sources is minimized. 

The existence of an optimal softening length follows from the assumed smoothness of 
the underlying mass distribution; no such optimal degree of smoothing can be defined for a 
"coUisional" code, where the accuracy of the force calculations can always be increased by 
decreasing the softening length. 

Our analysis will focus on the simplest A'"-body algorithm, a direct-summation code in 
which the smoothing is imposed via a fixed softening length associated with each particle. 
But the general principles carry over without significant change to codes that implement 
the smoothing in very different ways, e.g. via a grid, or via expansion of the potential in a 
truncated basis set. 



2. Definitions 

The gravitational force on particle i, 1 < i < N in a, direct-summation A'"-body code is 
given by 

where e, the softening length, determines the degree of smoothing. If the N particles 
represent an underlying, smooth distribution of mass, there also exists a "true" value 
Ffruei^i) for the force on particle i. For instance, if the N positions are generated from 
a Plummer density profile, Ffme is the gravitational force corresponding to the Plummer 
mass distribution. 

The optimal choice of e may be defined as the value that minimizes the average 
deviation between Fj and Ftmei'^i)- There are many possible ways to define this average 
deviation; the simplest is the average square error, or 

1 N 

ASE=-J2\Fi-Ftrue{^)\^. (2) 

i=l 
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We may treat F as a continuous function, defined via Eq. (1) at any point x. The 
continuous analog of the average square error is the integrated square error, or 

ISE = J p(x)|F(x) - Fi,,e(x)|'cix (3) 

with p(x) the true density, normahzed to unit total mass. 

If we imagine generating many A'"-particle Monte-Carlo realizations of the same smooth 
model, we can define the expectation or mean value of the ISE as 

MISE = E(ISE) = eJ p(x)|F(x) - F(x)t,,e|'c^x, (4) 

where E indicates an average over many reahzations. It is easy to show that the MISE 

contains contributions from two terms, representing the bias and the variance of the force 
estimates. The bias is the mean deviation of the computed force at some point from the 
true force, EF — Ffrue', its integrated square value is 

ISB = J p(x)|£;F(x) - Ft,„e(x)|'dx. (5) 

The variance is the mean square deviation of the force estimate from its mean value £^F; 
the integrated variance is IV = MISE — ISB. 

The bias increases with the softening length while the variance falls off - a greater 
amount of smoothing produces smaller average fluctuations in the forces but also tends 
to smooth over real features of small scale. For iV-particle realizations of a given model, 
therefore, there is an optimal choice of softening length Sopt such that the sum ISB + IV 
is minimized. One typically finds in problems of this sort (e.g. Scott 1992, p. 131) that 
the contribution to the total error from the bias and from the variance are of the same 
order when the smoothing length is optimized. Furthermore, the variance decreases with 
particle number while the bias is independent of N (e.g. Silverman 1986, p. 39); thus as N 
increases, one expects Cgpt to decrease. 



3. Two Examples 

We illustrate these ideas using N equal-mass points generated from two simple models. 
The first is a spherical Plummer model of unit scale length: 

The expectation value of the forces generated by a finite set of particles drawn from the 
density (6) will be spherically symmetric. We may therefore compute the terms in Eqs. (4) 
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and (5) on a linear grid in radius, taking averages over a large number of independent 3-D 
realizations of the particle distribution for any specified N. The number of Monte-Carlo 
reahzations was taken to be Nmc = 3 x e.g. 1000 reahzations for N — 3000. 

Fig. 1 shows how the contribution to the force errors from the bias and variance depend 
on € ior N — 1000. As expected, the bias increases with e while the variance decreases, 
producing an optimal value of e for which the MISE is minimized, Sopt ~ 0.16. At e = Sopt, 
roughly 60% of the total error comes from the variance and the remainder from the bias. 

Fig. 2 shows how the MISE vs. e curve changes with N. As N increases at fixed e, 
the variance decreases, roughly as A^~^. The result is an optimal e that decreases with N. 
Fig. 3 presents eopt and MISE(eopt) as functions of N. The two quantities vary roughly as 
power-laws in A^, with 

eopt ~ 1-1 X A^-°•2^ MlSEopt ~ 0.21 x AT-oee. (7) 



Our second example is based on a spherical model with the density law 

p(r) = i-r-^(l + r)-^ (8) 

which is a good approximation to the light distribution in early-type galaxies (Hernquist 
1990). Figure 3 shows Cppt and MISE(eopt) as functions of A^. The relations in this case can 
be approximated as 

eopt ~ 1-5 X MISEopt ~ 0.39 x N'^-^^ (9) 

The variation of topt with A^ in both models is such that the number of particles within one 
softening volume is approximately independent of A^. However in both cases the dependence 
seems to be significantly different than N~^^^. 



4. Discussion 

Selection of the softening length in A^-body codes is often made on the basis of the 
error due either to the bias or the variance of the force estimates. For instance, e may be 
chosen to be small enough that some feature like the core is resolved (bias), or large enough 
that two-body relaxation is not a problem (variance). Our discussion shows that the total 
error in the estimated forces in a coUisionless code is the sum of these two components, and 
that their contributions vary in opposite ways with e, leading to a formally optimal choice 
of e for any A". 
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The softening length in Eq. (1) may be interpreted as the window width of a kernel 
that is convolved with the density before computing the gravitational forces. Smoothing 
in coUisionless A^-body codes is usually carried out in other ways than via a kernel, e.g. 
on a grid. (Notable exceptions are the A^-body studies of galaxies based on the GRAPE 
hardware, which incorporates Eq. (1) [Ebisuzaki et al. 1993].) However, any smoothing 
algorithm may be interpreted in terms of a generalized kernel whose shape and width varies 
with position (Scott 1992, p. 155); in this sense, the analysis given here is general. One 
could use similar arguments to compute the optimal cell size in a grid-based algorithm, or 
the optimal number of terms in a basis-function expansion. 

It is well known that "coUisionless" A^-body codes often exhibit nearly as much 
relaxation as "coUisional" codes (e.g. Hernquist & Barnes 1990). The reason is that the 
variance in the force estimates - which, crudely speaking, is responsible for the relaxation - 
cannot be arbitrarily reduced without introducing a bias, as shown by Fig. 1. The primary 
design goal for coUisionless A'"-body experimenters (aside from speed) should be to smooth 
in such a way as to minimize the ISE of the force estimates at every time step, since by 
doing so they will have effectively minimized the variance. 

One way to decrease the ISE is to vary the softening length with position (and hence 
with time). The optimal way of doing this is well known when the quantity to be estimated 
is the density itself (Abramson 1982), but apparently nothing is known about the best way 
to vary e with position when estimating the gravitational forces. This would be a fruitful 
topic for further study. 

One would like to be able to estimate the optimal smoothing parameter directly from 
the A^ positions themselves, without having to know Ffrue- In a direct-summation code, 
automatic choice of e would presumably have to be based on a time-intensive bootstrap 
algorithm. More rapid choice of the smoothing parameter could be made in a code that 
represents the forces via sums of terms like Rim{T)Yim{0 , 4>) , with Yim an angular harmonic. 
The radial functions Rim may themselves be represented via basis-function expansions (e.g. 
Clutton-Brock 1973), but it seems to have escaped general notice that the job may be done 
more simply and efficiently using smoothing sphnes (Wahba 1990). Standard, and efficient, 
routines for computing the optimal degree of smoothing for such splines are widely available 
(e.g. Green & Silverman 1994). The result would be a code in which the radial dependence 
of the forces was computed with minimum ISE at every time step. 

One sometimes sees implementations in which N has been increased to very large 
values without corresponding adjustments in the degree of smoothing (e.g. Hernquist, 
Sigurdsson & Bryan 1995). Figure 2 shows that such a practice is likely to put the user far 
to the right of the minimum in the MISE curve. 
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Finally, we note that objective criteria for comparing the performance of different 
A^-body codes are generally lacking. A natural criterion for coUisionless codes would be the 
rate of decrease of the MISE with N. An optimal potential solver could be defined as one 
for which the MISE decreased most rapidly with A^. 

This work was supported by NSF grant AST 90-16515 and NASA grant NAG 5-2803. 
I thank G. Quinlan for helpful discussions. 
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Fig. 1. — Bias-variance tradeoff in the force evaluation errors for N — 1000 realizations of a 
spherical Plummer model. 




Fig. 2.— MISE(e) for N = 30, 100, 300, 1000, 3000, 10000 and 30000; N increases downward. 
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Fig. 3.- 

Dependence of the optimal softening length €opt, and the MISE at e = e^pt, on number of 
particles N in two models: Eq. (6) (squares) and Eq. (8) (circles). The relations of Eqs. 

(7) and (9) are indicated by solid lines. 



